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Abstract We propose two algorithms to provide a full preliminary orbit of an Earth- 
orbiting object with a number of observations lower than the classical methods, such 
as those by Laplace and Gauss. The first one is the Virtual debris algorithm, based 
upon the admissible region, that is the set of the unknown quantities corresponding to 
possible orbits for objects in Earth orbit (as opposed to both interplanetary orbits and 
ballistic ones) . A similar method has already been successfully used in recent years for 
the asteroidal case. The second algorithm uses the integrals of the geocentric 2-body 
motion, which must have the same values at the times of the different observations for 
a common orbit to exist. We also discuss how to account for the perturbations of the 
2-body motion, e.g., the J2 effect. 

Keywords Space debris ■ Orbit determination • Admissible region • Keplerian 
integrals 



1 Introduction 

The near-Earth space, filled by more than 300000 artificial debris particles with diam- 
eter larger than 1 cm, can be divided into three main regions: the Low Earth Orbit 
(LEO), below about 2000 km, the Medium Earth Orbit (MEO), above 2000 km and 
below 36000 km, and the Geosynchronous Earth Orbit (GEO) at about 36000 km of 
altitude. Currently the orbits of more than 12000 objects larger than about 10 cm are 
listed in the so called Two Line Elements (TLE) catalogue. To produce and maintain 
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such a catalogue a large number of optical and radar observations are routinely per- 
formed by the United States Space Surveillance Network. Nowadays also Europe has 
launched its Space Situational Awareness (SSA) initiative aimed to increase the knowl- 
edge of the circumterrestrial environment. In this context the availability of efficient 
methods and algorithms for accurate orbit determination is extremely important. 

Given two or more sets of observations, the main problem is how to identify which 
separate sets of data belong to the same physical object (the so-called correlation 
problem) . Thus the orbit determination problem needs to be solved in two stages: first 
different sets of observations need to be correlated, then an orbit can be determined; 
this combined procedure is called linkage in the literature Milani 1999 . 

In this paper we describe two different linkage methods, for both optical and radar 
data. By using the attributable vector (Sec. [2} we summarize the information contained 
in either optical or radar data. In Sec. [3] we describe the admissible region and the Vir- 
tual debris algorithm Tommci ct al. 2007] and we propose a general scheme to classify 
observed objects. Sec. [4] deals with the Keplerian integrals method, first introduced by 
Gronchi ct al. 2009 for the problem of asteroid orbit determination. Furthermore, the 
inclusion of the effect due to the non-spherical shape of the Earth is discussed. Finally, 
in Sec. [5] a sketch of the general procedure for the full process of correlation of different 
observations is outlined. 

2 Observations and attributables 

Objects in LEO are mostly observed by radar while for MEOs and GEOs optical sensors 
are used. In both cases, the batches of observations which can be immediately assigned 
to a single object give us a set of data that can be summarized in an attributable, that 
is a 4-dimensional vector. To compute a full orbit, formed by 6 parameters, we need to 
know 2 further quantities. 

Thus the question is the identification problem, also called correlation in the debris 
context: given 2 attributables at different times, can they belong to the same orbiting 
object? And if this is the case, can we find an orbit fitting both data sets? 

Let (p, a, S) G K + x [0, 27r) x (— ty/2, n/2) be spherical coordinates for the topocentric 
position of an Earth satellite. The angular coordinates (a, 8) are defined by a topocen- 
tric reference system that can be arbitrarily selected. Usually, in the applications, a is 
the right ascension and S the declination with respect to an equatorial reference system 
(e.g., J2000). The values of range p and range rate p are not measured. 

We shall call optical attributable a vector 

Aopt = {a, S, a, 6) € [0, 2tt) x (tt/2, tt/2) x R 2 , 

representing the angular position and velocity of the body at a time t in the selected 
reference frame. 

Active artificial satellites and space debris can also be observed by radar: because 
of the 1/p dependence of the signal to noise for radar observations, range and range- 
rate are currently measured only for debris in LEO. When a return signal is acquired, 
the antenna pointing angles are also available. Given the capability of modern radars 
to scan very rapidly the entire visible sky, radar can be used to discover all the debris 
above a minimum size while they are visible from an antenna, or a system of antennas. 

When a radar observation is performed we assume that the measured quantities (all 
with their own uncertainty) are the range, the range rate, and also the antenna pointing 
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direction, that is the debris apparent position on the celestial sphere, expressed by two 
angular coordinates such as right ascension a and declination S. The time derivatives 
of these angular coordinates, a and 8, are not measured. 
We define radar attributable a vector 

A rad = (a,5,p,p) 6 [-tt.tt) x (-jr/2, jt/2) x E+ x K, 

containing the information from a radar observation, at the receive time t. 

To define an orbit given the attributable A we need to find the values of two 
unknowns quantities (e.g., p and p in the optical case, a and S in the radar case), that, 
together with the attributable, give us a set of attributable orbital elements: 

X = [a,5,d,8,p,p] 

at a time i, computed from t taking into account the light-time correction: t = t — p/c. 
The Cartesian geocentric position and velocity (r, r) can be obtained, given the observer 
geocentric position q at time t, by using the unit vector p = (cos a cos 8, sin a cos 5, sin 8) 
in the direction of the observation: 

dp dp . „ • 
r = q + , r = q + pp + p— , — = ap a + 8p s , 

p a — (— sin a cos S, cos a cos 5, 0) , Ps = (— cos a sin 5, — sin a sin 5, cos 8) . 



3 Admissible region theory 

Starting from an attributable, we would like to extract sufficient information from it 
in order to compute preliminary orbits: we shall use the admissible region tool, as 
described in Tommci ct al. 2007 . For easy of reading we recall here the basic steps of 
the theory. 

The admissible region replaces the conventional confidence region as defined in the 
classical orbit determination procedure. The main requirement is that the geocentric 
energy of the object is negative, that is the object is a satellite of the Earth. 



3.1 Optical admissible region 

Given the geocentric position r of the debris, the geocentric position q of the observer, 
and the topocentric position p of the debris we have r = p + q. The energy (per unit 
of mass) is given by 



2" ™ ||r(p)||' 

where p is the Earth gravitational parameter. Then a definition of admissible region 
such that only satellites of the Earth are allowed includes the condition 

£(P,P)<0 (1) 

that could be rewritten as 

2£(p, p)=p 2 + c lP + T(p) - -^§= < , 
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T(p) = c 2 p 2 + c 3 p + c 4 , S(p) = p 2 + c 5 p + c 
and coefficients Cj depending on the attributable ITommei et al. 2007] : 

eo = ||q|| 2 , ei=2q-/5, ci = a 2 cos 2 <5 + 8 2 — rj 2 , 

C3 = 2(dq- p tt + <5q- p s ) , c 4 = ||q|| 2 , c 5 = 2 q ■ p , 

where r\ is the proper motion. In order to obtain real solutions for p the discriminant 
of 2£ (polynomial of degree 2 in p) must be non-negative: 

Z\=^i-T(p) + -^>0. 

This observation results in the following condition on p: 

-^=>Q(p)=c 2 p 2 +c 3 p + 1 , 1 = c 4 ~^. (2) 

Condition ((2} can be seen as an inequality involving a polynomial V(p) of degree 6: 

V(p) = Q 2 (p)S(p)<4p 2 . 

Studying the polynomial V(p) and its roots, as done by [Milani et al. 2004] . the con- 
clusion is that the region of (p, p) such that condition fl]) is satisfied can admit more 
than one connected component, but it has at most two. In any case, in a large number 
of numerical experiments with objects in Earth orbit, we have not found examples with 
two connected components. 

The admissible region needs to be compact in order to have the possibility to 
sample it with a finite number of points, thus a condition defining an inner boundary 
needs to be added. The choice for the inner boundary depends upon the specific orbit 
determination task: a simple method is to add constraints p m i n < p < pmax allowing, 
e.g., to focus the search of identifications to one of the three classes LEO, MEO and 
GEO. Another natural choice for the inner boundary is to take p > h a tm where h a t m 
is the thickness of a portion of the Earth atmosphere in which a satellite cannot remain 
in orbit for a significant time span. As an alternative, it is possible to constrain the 
semimajor axis to be larger than R§ + h a tm = r m i n , and this leads to the inequality 

£{P,p) > -r— = £rnin , (3) 

which defines another degree six inequality with the same coefficients but for a different 
constant term. The qualitative structure of the admissible region is shown in Fig. \T\ 

Another possible way to find an inner boundary is to exclude trajectories impacting 
the Earth in less than one revolution, that is to use an inequality on the perigee rp, 
already proposed in [Maruskin et al. 2009 : 

r P = a(l - e) > r mm . (4) 

Note that this condition naturally implies (|3]) and p > h a t m - To analytically develop 
the inequality Q we use the 2-body formulae involving the angular momentum: 

c(p,p) = r x r = Dp + Ep 2 + Fp + G (5) 
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D = q x p, E = p x (dp a + Sp s ), 

F = q X (ap a + Sp s ) + q X p, G = q X q 

and substituting in Q we obtain: 



1 + 2g I^H 2 < 1 + 2£> """ (g) 
^ /i 

Since the left hand side is e > 0, we need to impose 1 + 2£r m i n /p > 0: this is again 
a > r min- By squaring ((SJl we obtain: 

|c| | > 27" m 2n(/^ ~t~ &Tmin) ■ 

The above condition is an algebraic inequality in the variables (p, p): 

(r^„ - I |D|| 2 ) P 2 - P(p)p - U( P ) + r 2 mm T{p) - ^kp^ < , (7) 



P(p) = 2D ■ Ep 2 + 2D Fp + 2D G , 

U(p) = ||E|| 2 p 4 + 2E ■ Fp 3 + (2E ■ G + ||F|| 2 )p 2 + 2F ■ Gp + ||G|| 2 - 2r min p. 

The coefficient of p 2 is positive, thus to obtain real solutions for p the discriminant of 
(0 must be non negative: 

A F = P 2 (p) + 4(r 2 nm - j|D|| 2 ) (u(p) + T 2 mm T{p) + > • 
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This condition is equivalent to the following: 

2 M > w{p) _ < r Ln ~ I|D|| 2 )(^(P) + r 2 mm T(p)) + P 2 (p) ^ (g) 



Note that the inequality (|8} is similar to ([2]). However, in this case, the function in 
the right hand side is much more complicated, and there is no easy way to use the 
condition Q to explicitly describe the boundary of the admissible region; e.g., we do 
not have a rigorous bound on the number of connected components. This condition Q 
will be used only a posteriori as a filter (Sec. 13 .3 P . 




Fig. [2] shows also this inner boundary; note that the boundaries of the regions 
defined by ([3]) and by p > h a t m are also plotted in the figure, but these constraints 
are not necessary. We have also plotted an alternative outer boundary constraining the 
apocenter r ^ at some large value r m ax- 



) ^rmax 
I 1 1 ^1 1 ^ Qr-max ~t~ ^Tmax) 

this outer boundary can be used in the same way, as an a posteriori filter. 



3.2 Radar admissible region 

Given a radar attributable A ra d, we define as radar admissible region for a space debris 
the set of values of (a, 8) such that 



£ (a, 5) — z\a 2 + Z2$ 2 + z%a + Z46 + Z5 < , 



(9) 
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E=0 




da/dt cos5 (rad/day) 

Fig. 3 An example of admissible region, radar case, in the (a cos 8,5) plane. The region 
(painted in grey) is the circular annulus bounded by the two level curves of the energy [E = 
E m in) and (E = 0). 



where 2y depend on the attributable Tommci et al. 2007 : 



2 2 c 

z\ = p cos o , 
z 4 = pq- ps/2, 



Z2 = p , 

.2 . 2u 

Z5 — P + Cip + c 4 ■ 



Z3 = PQ ■ Pa/2 , 



The boundary of the admissible region is then given by £(d, 5) = and this equation 
represents an ellipse with its axes aligned with the coordinate axes in the (a, 6) plane. 
Actually, in a plane (d cos S, 5), with the axes rescaled according to the metric of the 
tangent plane to the celestial sphere, the curves £ (a, S) = constant are circles. 

The region defined by negative geocentric energy, the inside of a circle, is a compact 
set, and the problem of defining an inner boundary is much less important than in 
the optical attributable case. Anyway, it is possible to define an inner boundary by 
constraining the semimajor axis a > r m j„, that is by eq. (J3J), resulting in a concentric 
inner circle, thus in an admissible region forming a circular annulus (see Fig. [5J • 

It is also possible to exclude the ballistic trajectories by imposing the condition Q 
in which d, S are to be considered as variables. The angular momentum is given by 



c = r x r = Ad + B<5 + C , 



(10) 



A = pr x p a , B = p r X pg , C = rxq + pqX/3. 
The condition on the pericenter is expressed by a polynomial inequality of degree 2: 

w\a 2 + W2a5 + W38 2 + w^a + W56 + w q > , 



W2 



2A B, 



w 3 
W4 



l|B|. 
2(A- C- 



2rt 



w 5 = 2(B ■ C - r,, 



, z i) 



w 6 



2r 
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Fig. 4 An example of admissible region, with the further condition on the pericenter distance 
(rp > r m j n ), bounded by an ellipse. 



150 - 




-150 ^ , , , , , , , 

-150 -100 -50 50 100 150 

da/dt cos5 (rad/day) 

Fig. 5 An example of admissible region, with the further condition on the pericenter distance 
(rp > r m j n ), bounded by an hyperbola. 



Thus the admissible region can be geometrically described as a region bounded by 
three conies: the first two are concentric circles, the third one can be either an ellipse 
or an hyperbola (depending on the sign of wiw^ — «;|/4), with a different center and 
different symmetry axes. Fig. [4] and [5] show the possible qualitatively different cases. 
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Fig. 6 An example of admissible region, defined by imposing negative geocentric energy, 
for an optical attributable, with the Dclaunay triangulation. The nodes of the triangulation 
corresponding to the ballistic trajectories (on the left of the curve cutting the outer part of 
the triangulation) can be discarded. 




doVdt cos8 



Fig. 7 An example of admissible region, defined by imposing £ m i n < £ < 0, for a radar 
attributable, with the cobweb sampling. The nodes of the cobweb corresponding to the ballistic 
trajectories (between the two branches of the hyperbola) can be discarded. 
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3.3 Virtual debris algorithm 

The admissible region can be used to generate a swarm of virtual debris: we sample 
it using the Delaunay triangulation Mil ani et al. 2004] for the optical case and the 
cobweb Tommei et al. 2007] for the radar shown in Fig. [6] and [7] The condi- 

tion on the pericenter is not used at this step, because we could lose some important 
geometrical properties: this condition is used as filter, the nodes with a low pericenter 
are discarded. 

The idea is to generate a swarm of virtual debris Xi, corresponding to the nodes 
of the admissible region of one of the two attributables, let us say Ai . Then we com- 
pute, from each of the Xi, a prediction Ai for the epoch ti, each with its covariance 
matrix Fj^. . Thus for each virtual debris Xi we can compute an attribution penalty 
K\ [Milani and Gronchi 2009] and use the values as a criterion to select some of the 
virtual debris to proceed to the orbit computation. 

Thus the procedure is as follows: we select some maximum value K m ax for the 
attribution penalty and if there are some nodes such that K\ < K ma x we proceed to 
the correlation confirmation. If this is not the case, we can try with another method, 
such as the one described in Sec. H] 



3.4 Universal classification of objects 




Fig. 8 Partitioning of the (p, p) half plane p > in regions corresponding to different popula- 
tions, for an optical attributable with proper motion r] = 10.1980. The labels mean: L Launch, 
R Reentry, ES Earth Satellite, A Asteroid, ISC Interstellar Comet, ETA ET Arriving, ETL 
ET Leaving, IPR Interplanetary Reentry, IPL Interplanetary Launch, rad/day. 



11 



400 



ISC 




-600 



-600 



-400 



-200 





da/dt cos5 (rad/day) 



200 



400 



600 



Fig. 9 Partitioning of the (a cos 5,(5) plane in regions corresponding to different populations, 
for a radar attributable with p = 1 R^. The labels mean: L/R Launch or Reentry, ES Earth 
Satellite, A Asteroid, ISC Interstellar Comet, ET interstellar launch/reentry, IPL/IPR Inter- 
planetary Launch or Reentry. 

The method of the admissible region is also useful to provide insight on the re- 
lationship between the different populations, in particular how they can mix in the 
observations. For a given optical attributable, supposedly computed from a short arc 
of optical observations, the Fig. [8] shows the region in the (p, p) half-plane p > where 
Earth satellites (ES) can be, but also where ballistic trajectories (either launches L or 
reentries R) can be, and where an asteroid serendipitously found in the same obser- 
vations would be. Other more exotic populations, which are very unlikely, also have 
their region in the half plane: e.g., there are regions for direct departure/arrival to the 
Earth from interstellar space, which we have labeled as ET trajectories. 

The same "universal" figure can be generated from a given radar attributable 
(Fig [9]). In this case the regions corresponding to different populations partition the 
plane (a cos S, 5). The curve £ S un = 0, for the heliocentric energy, has been computed 
with formulas very similar to the ones for the geocentric energy. 

4 Keplerian integrals method 

We shall describe a method proposed for the asteroid case in [Gronchi et al. 2009| 
and based on the two-body integrals, to produce preliminary orbits starting from two 
attributables A\ , A2 of the same object at two epochs t\ , ti. We assume that the orbit 
between t\ and t<x is well approximated by a Keplerian 2-body orbit, with constant 
energy £ and angular momentum vector c: 



{ 



B(ti) - £(ta) = 
c(ti) - c(t 2 ) = 



(11) 
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4.1 Optical case 

Using ([5]| and by scalar product between with the first equation of (|11[1 and Di x D2 
we obtain the scalar equation of degree 2: 

(Di x D 2 ) ■ (ci - c 2 ) = q(pi, P2)=0- 

Geometrically, this equation defines a conic section in the (pi, p 2 ) plane. By the formu- 
lae giving pi, p2 as a function of p\, p2 derived from the angular momentum equations: 

. _ (E 2 p| + F2P2 + G 2 - Eipf - Fipi - Gi) x D 2 
Pl HDxxDall 2 
. _ Di x (Eip| + Fipi + Gj_- E 2P | - F 2 p 2 - G 2 ) 
92 ||DixD 2 || 2 

the energies £\, £2 can be considered as functions of pi, p 2 only. Thus we obtain: 

£i(pi,P2) - £-i{p\,pi) = 
q(pi,P2) = 

a system of 2 equations in 2 unknowns, already present in [Taff and Hall 1977] : they 
proposed a Newton- Raphson method to solve the system, but this results into a loss of 
control on the number of alternate solutions. In Gronchi ct al. 2009 the authors have 
applied the same equations to the asteroid problem, and proposed a different approach 
to the solution of the system. 

The energy equation is algebraic, but not polynomial, because there are denomi- 
nators containing square roots. By squaring twice it is possible to obtain a polynomial 
equation p(pi, p 2 ) = 0: the degree of this equation is 24. Thus the system 

P{P1,P2) = 

q(pi,P2) = 

has exactly 48 solutions in the complex domain, counting them with multiplicity. Of 
course we are interested only in solutions with pi, p 2 real and positive, moreover the 
squaring of the equations introduces spurious solutions. Nevertheless, we have found 
examples with up to 11 non spurious solutions. 

We need a global solution of the algebraic system of overall degree 48, providing at 
once all the possible couples (pi,p 2 ). This is a classical problem of algebraic geometry, 
which can be solved with the resultant method: we can build an auxiliary Sylvester 
matrix, in this case 22 x 22, with coefficients polynomials in p 2 , and its determinant, 
the resultant, is a polynomial of degree 48 in p 2 only. The values of p 2 appearing in 
the solutions of the polynomial system are the roots of the resultant [Cox et al. 1 996 . 

The computation of the resultant is numerically unstable, because the coefficients 
have a wide range of orders of magnitude: we have to use quadruple precision. Once 
the resultant is available, there are methods to solve the univariate polynomial equa- 
tions, providing at once all the complex roots with rigorous error bounds [Bini 1996| . 
Given all the roots which could be real, we solve for the other variable pi, select the 
positive couples (pi,p 2 ) and remove the spurious ones due to squaring. If the number 
of remaining solutions is 0, the attributables cannot be correlated with this method. 
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4.2 Radar case 

The formulae for geocentric energy and angular momentum are given by (|9} and (|10[) , 
polynomials of degree 2 and 1 in the unknowns (a, 5), respectively. The system 
has overall algebraic degree 2: such a system can be solved by elementary algebra. 
The angular momentum equations are 

Aidi + Bi<5i + Ci = A 2 A 2 + B 2 <5 2 + C 2 (12) 

that is a system of 3 linear equations in 4 unknowns (ax, 5\, d 2 , <5 2 ) and can be solved 
for three unknowns as a function of one of the four. For example, by scalar product 
between (|12l) and Bi x A 2 we have 

. _ A 2 ■ (Bi x B 2 )S 2 - (Ci - C 2 ) ■ (A 2 x Bi) 
ai Bi • (Ai x A 2 ) 

and in a similar way we obtain 

■ _ B 2 ■ (Ai x A 2 ) 6 2 - (Ci - C 2 ) ■ (Ai x A 2 ) 
1 Bi ■ (Ai x A 2 ) 

. _ Ai ■ (Bi x B 2 )<5 2 - (Ci - C 2 ) ■ (Ai x BQ 
Q2 Bi • (Ai x A 2 ) 

When the equations for, say, (di,d 2 ,c5i) as a function of <5 2 are substituted in 
the equation for the energies £i(di,<5i) = £ 2 (d 2 ,<5 2 ) we obtain a quadratic equation 
in <5 2 , which can be solved by elementary algebra, giving at most two real solutions. 
Geometrically, equation (|12[) can be described by a straight line in a plane, e.g., in 
(d 2 , <5 2 ), where the energy equation defines a conic section. 



4.3 Singularities 

There are some cases in which the Keplerian integrals method can not be applied. 

In the optical case we have to avoid the condition Di x D 2 = (qi x pi) x (q 2 x p 2 ) = 
0. This can happen when: 

— qi is parallel to p\, i.e. the observation at time t\ is done at the observer zenith; 

— q 2 is parallel to p 2 , i.e., the observation at time t 2 is done at the observer zenith; 

— qi, q 2 , pi and p 2 are coplanar. This case arises whenever a geostationary object 
is observed from the same station at the same hour of distinct nights. 

As it is normal, the mathematical singularity is surrounded by a neighborhood in which 
the method is possible for zero error (both zero observational error and zero rounding off 
in the computation), but is not applicable in practice. E.g., for nearly geosynchronous 
orbits, even if they are not geostationary, and for hours of observations in different 
nights different by few minutes, this method fails. 

In the radar case the procedure fails only if the four vectors Ai, A 2 , Bi and B 2 
do not generate a linear space of dimension 3, i.e., when: 

Ar(BixB 2 )=0 
Br(AixA 2 ) = 1 ' 
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For i = 1 we obtain 

\plP2[P52 ■ (ri x r 2 )][ri ■ (p a l X Psi)] = 
\p?P2[/3q2 • (ri x r 2 )][ri ■ (p si x p al )] = 

and for i — 2 the formulae are analogous. Thus there is singularity when: 

— ri is parallel to r 2 ; 

— r-j ■ {psi x p Q j) = cos5j(qj ■ pi + pi) = 0, but this can never happen, apart from 
coordinate singularities, because <li ■ Pi > 0; 

— r l ■ {Psi x Pai) = and n • (p si x p Si ) = 0, i.e., p m and p Si for i = 1, 2 belong to 
the orbital plane. 



4.4 Preliminary orbits 

Once a solution of JTT} is computed the values of attributable elements can be obtained 
for the epochs t\ and ii, and they can be converted into the usual Keplerian elements: 

(a j , ej , Ij , Hj , coj ,lj) , j = 1, 2 , 

where ^ are the mean anomalies. The first four Keplerian elements (a,j, ej, Ij, Qj) are 
functions of the 2-body energy and angular momentum vectors Ej, Cj, and are the 
same for j = 1, 2. Thus the result can be assembled in the 8-dimensional vector: 

H=(V,* 1 ,* 2 ) , V = (o,e,/,fl), * 1 = (wi,/i), #2 = (wa,/2). (13) 

There are compatibility conditions between <fr\ and ^2 to be satisfied if the two at- 
tributables belong to the same object: 

wi =w 2 , ^1 =l2+n(tx-h) , (14) 

where n = n(a) is the mean motion. We cannot demand the exact equality in the 
formulae above, because of various error sources, including the uncertainty of the at- 
tributable, and the changes on the Keplerian integrals due to the perturbations with 
respect to the 2-body model. Thus we need a metric to measure in an objective way 
the residuals in the compatibility conditions. 



4.5 Covariance propagation 

The two attributables Ai , Ai used to compute the coefficients of equations have 
been computed from the observations by using a least squares fit to the individual 
observations, thus 4x4 covariance matrices and r^ 2 are available; they can be 
used to form the block diagonal 8x8 covariance matrix for both attributables r^. 
The Keplerian integral method allows to compute explicitly the vector H of (|13|) and, 
by means of the implicit function theorem, its partial derivatives, thus it is possible 
by the standard covariance propagation formula Milani and Gronchi 2009 [Sec. 5.5] to 
compute also /jj, the covariance of H. With another transformation we can compute 
the average elements <I>q = (<Pi + $2)/2 (as the best value for the angular elements at 
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time to = (ti + t2)/2) and the discrepancy A$ in the compatibility conditions (|14|) . 
and to propagate the covariance also to this 8-dimensional vector: 

Fa Fh rv,<p ,A& ■ 

The above argument is a generalization of the one in [Gronchi et al. 2009] . where ex- 
plicit computations are given for the optical attributables case. 

In the 8x8 covariance matrix -TV,* ,4#' * ne l° wer right 2x2 block is the marginal 
covariance matrix of A<P, from which we can compute the normal matrix and the \ 2 '- 

which can be used as control, that is the discrepancy in the compatibility conditions is 
consistent with the observation error and the correlation between the two attributables 
is considered possible only if Xa& — Xmax- 

The upper left 6x6 block is the covariance matrix of the preliminary orbit, that is 
of the orbital elements set (V, <&o) (at epoch to). Although this preliminary orbit is just 
a 2-body solution, it has an uncertainty estimate, arising from the (supposedly known) 
statistical properties of the observational errors. This estimate neglects the influence of 
perturbations, such as the spherical harmonics of the Earth gravity field, the lunisolar 
differential attraction and the non-gravitational perturbations; nevertheless, if the time 
span I2 — i\ is short, the covariance obtained above can be a useful approximation. 



4.6 Precession model 



We can generalize the method, including the effect due to the non-spherical shape of 
the Earth. The averaged equation for Delaunay's variables £, g = w, z = Q, L — ^/Jm, 
G = Wl - e 2 and Z = Gcos/ are |Roy 2005| [Sec. 10.4]: 



-JV 



Ra 



q = —n 
y 4 



L = G 



a 

a 

Z = 



J 2 (l -3 cos 2 I) 

(1-e 2 ) 3 / 2 

J 2 (4- 5 sin 2 J) 

(1-e 2 ) 2 

2 t r 
J2 COS I 



(15) 



(1-e 2 ) 



2v2 



where J 2 is the coefficient of the second zonal spherical harmonic of the Earth gravity 
field. To apply in this case the Keplerian integrals method, we can not use the equations 
assuming conservation of the angular momentum. From (|15[) we can replace (|11|) with: 



£i = £a 

Cl • Z = C 2 • Z 

II Il2 II |2 
ll c l|| = ll c 2|| 

cos(«2) = cos(zi) cos(i(?2 — ii)) 



(16) 



i(ai)sin(i(t2 - *i)) 



In the optical case the first equation is algebraic and by squaring twice is possible to 
obtain a polynomial equation; in the radar case this relation is already polynomial. The 
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second and the third equations are always polynomial, while the last equation needs 
to be linearized in the parameter z(ti — t%): 

cos 22 = cos z\ — i(t2 — t\) sin z\ . (17) 

The following relationships hold: 

ZXCj . . ZXCj , . 2\ !|ci|| 2 
COS Zj = jp rr ■ X , S1I1Z, = yp; rr • y , Cl(l — 6 J = , 



cos/=^— ^, n=\l-— , &(t 2 -ti) 

\\ c i\\ V A* 



8£f — £*/-8£? (ci • *) 

where £ = —Sfi^R^fo ~ Substituting in (JXTJ we obtain 



ZXC2 , ZXCl . ''V "'"'I v 1 I Z X Cl 



ZJ-8£f{a -z) 



| |z X C2 1 1 ||z X Cl|| ll c l|| 5 |z X Cl 1 

Since llz x ell = II ell sin / is constant we have: 



|ci|| 5 [z x (c 2 - ci)] ■ x = -£yJ-8£% (ci ■ z)(z x ci) • y , 

that is an algebraic equation. Furthermore, by squaring twice in the optical case and 
only once in the radar case it is possible to obtain a polynomial equation. 

Finally the new compatibility conditions, in place of (|14[) need to take into account 
the precession of the perigee and the secular perturbation in mean anomaly: 

.91 = .92 + Kh - h) , h=h+ - h) ■ 

The overall degree of system (|16[1 is summarized in Table [T] We conclude that this 
method is unpractical for optical attributables, could be used for radar attributables, 
with computational difficulties comparable with the optical case without precession. 



Table 1 Degrees of the equations in system II16H 





Optical case 


Radar case 


Si =s 2 


16 


2 


Cl ■ Z = C2 ■ z 


2 


1 


l|ci|| 2 = ||c 2 || 2 


4 


2 


cos 22 = cos 21 i(<2 — ti) sin 2i 


54 


12 


Total 


6912 


48 



To solve the problem (even in the optical case) we begin by considering the para- 
metric problem i = K, where K is constant. Thus we replace (|lip with: 

£ i - £2 = 
Rd - R T c 2 = 

where R is the rotation by AS1/2 — K(t% —t\)/2 around z. This means that for a fixed 
value of K the problem has the same algebraic structure of the unperturbed one. The 
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only thing needed is to substitute Di, Ei, Fi and Gi with -RDi, REi, RFi and RGi 
in the optical case and Ai, Bi and Ci with RAi, RBi and RCi in the radar case; 
similarly the vectors with index 2 are multiplied by R T . 

The compatibility conditions contain the precession of the perigee and the secular 
perturbation in mean anomaly, related to the one of the node by linear equations 

91 = 92 + KC g (t! - h) , h=h + {n + KC t )(h - h) , 

where the coefficients Cg,Ci can be easily deduced from (|15[1 . Thus we can compute 
the Xa&(K) an d set up a simple procedure to minimize this by changing K, then the 
control on the acceptability of the preliminary orbit is 

minxLfW < Xmax ■ 

5 Correlation confirmation 

The multiple orbits obtained from the solutions of the algebraic problem are just pre- 
liminary orbits, solution of a 2-body approximation (as in the classical methods of 
Laplace and Gauss), or possibly of a J2-only problem. They have to be replaced by 
least squares orbits, with a dynamical model including all the relevant perturbations. 

Even after confirmation by least squares fit, it might still be the case that some 
linkages with just two attributables can be false, that is the two attributables might 
belong to different objects. This is confirmed by the tests with real data reported in 
Tommei et al. 2009] for the Virtual debris method and in [M ilani e t al. 2009] for the 
Keplerian integrals method. Gronchi ct al. 2009 have found the same phenomenon in 
a simulation of the application of the same algorithm to the asteroid case. Thus every 
linkage of two attributables needs to be confirmed by correlating a third attributable. 

The process of looking for a third attributable which can also be correlated to 
the other two is called attribution [Milani 1 999, Mila ni et al. 2001) . From the available 
2-attributable orbit with covariance we predict the attributable Ap at the time t% of 
the third attributable, and compare with ^3 computed from the third set of obser- 
vations. Both Ap and A3 come with a covariance matrix, we can compute the x 2 °f 
the difference and use it as a test. For the attributions passing this test we proceed 
to the differential corrections. The procedure is recursive, that is we can use the 3- 
attributable orbit to search for attribution of a fourth attributable, and so on. This 
generates a very large number of many-attributable orbits, but there are many dupli- 
cations, corresponding to adding them in a different order. 

By correlation management we mean a procedure to remove duplicates (e.g., A — 
B — C and A — C = B) and inferior correlations (e.g., A — B = C is superior to both 
A = B and to C — D, thus both are removed). The output catalog after this process 
is called normalized. In the process, we may try to merge two correlations with some 
attributables in common, by computing a common orbit Milani ct al. 2005 . 

6 Conclusions 

We have described two algorithms to solve the linkage problem, that is to compute 
an orbit for an Earth-orbiting object observed in two well separated arcs. The first 
method exploits the geometric structure of the admissible region of negative geocentric 
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energy orbits, which is sampled to generate virtual orbits. The latter are propagated 
in time to find other observations which could belong to the same object. The second 
method exploits the integrals of the 2-body problem, which are constant even over a 
significant time span and thus should apply to both observed arcs of the same object. 

This top level description is enough to understand that the Virtual debris algo- 
rithm should be applied to short time intervals between observed arcs, less than one 
orbital period or at most a few orbital periods. The Keplerian integrals method can 
be used for longer time spans, spanning several orbital periods; it is near to a sin- 
gularity for very short time spans and in some other near-resonance conditions, such 
as observations of a geosynchronous orbits at the same hour in different nights. We 
conclude that each method should be used in the cases in which it is most suitable. 
Both algorithms have been tested for the optical case with real data from the ESA Op- 
tical Ground Station Tommci et al. 2009 , Milani ct al. 2009 with good results. The 
analogous algorithms have been tested for asteroids in simulations of next generation 
surveys [Milani et al. 20 05 . Gronchi ct al. 2009 . Future work should include the tests 
of the radar case and the solution of other related problem, like orbit identification 
between two objects for which an orbit is already available. 
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